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Abstract 

A modification scheme to the ensemble Kalman filter (EnKF) is introduced based 
on the concept of the unscented transform (Julier et al., 2000; Julier and Uhlmann, 
2004), which therefore will be called the ensemble unscented Kalman filter (EnUKF) 
in this work. When the error distribution of the analysis is symmetric (not necessar- 
ily Gaussian), it can be shown that, compared to the ordinary EnKF, the EnUKF 
has more accurate estimations of the ensemble mean and covariance of the back- 
ground by examining the multidimensional Taylor series expansion term by term. 
This implies that, the EnUKF may have better performance in state estimation than 
the ordinary EnKF in the sense that the deviations from the true states are smaller. 
For verification, some numerical experiments are conducted on a 40-dimensional sys- 
tem due to Lorenz and Emanuel (Lorenz and Emanuel, 1998). Simulation results 
support our argument. 
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1 Introduction 



The Kalman filter (KF) is a recursive data processing algorithm [23]. It op- 
timally estimates the states of linear stochastic systems that are driven by 
Gaussian noise, and are observed through linear observation operators, which 
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possibly also suffer from additive Gaussian errors. However, if there exists non- 
linearity from either the dynamical systems or the observation operators, or, if 
neither the dynamical noise nor the observational noise follows any Gaussian 
distribution, then the Kalman filter becomes suboptimal. To tackle the prob- 
lems of nonlinearity and non-Gaussianity, there are some strategies one may 
employ. For example, to handle the problem of nonlinearity, one may expand 
the nonlinear function in a Taylor series locally and keep the expansion terms 
only up to second order. This leads to the extended Kalman filter (EKF) (e.g., 
[7]). To deal with the problem of non-Gaussianity, one may specify a Gaus- 
sian mixture model (GMM) to approximate the underlying probability density 
function (pdf), such that the KF algorithm is applicable to the individual dis- 
tributions of the GMM [27]. More generally, one may adopt the sequential 
Monte Carlo method (also known as the particle filter, e.g., |32j), which uti- 
lizes the empirical pdf obtained from a number of particles to represent the 
true pdf, wherein the problems of both nonlinearity and non-Gaussianity are 
taken into account during the pdf approximation. 

For practical large-scale problems like weather forecasting, the computational 
cost is another issue of great concern. In such circumstances, direct application 
of the KF or EKF scheme is prohibitive because of the computational cost of 
evolving the full covariance matrices forward. While for the particle filter, 
because of its slow convergence rate, the required number of the particles 
for proper approximations may be well above many thousands for even low 
dimensional nonlinear systems. 

For the sake of computational efficiency, the ensemble Kalman filter (EnKF) 
was proposed [8]. It is essentially a Monte Carlo implementation of the Kalman 
filter. At the beginning of each assimilation cycle, it is assumed that one 
has an ensemble of the system states, called the background ensemble, which 
can usually be obtained from the previous assimilation cycle. Then, given 
additional information from the observations, one applies the KF scheme to 
update each individual member of the background ensemble. To do this, the 
mean and error covariance of the background are approximated by the sample 
mean and covariance of the ensemble. After the updates, one will obtain a 
new set of the system states, called the analysis ensemble in this work, which 
will then be used to estimate the true mean and covariance of the underlying 
system state. By propagating the ensemble members of the analysis forward 
through the system model, one obtains a new ensemble of the background for 
the next assimilation cycle. Therefore, by using the ensembles to approximate 
the true statistics (e.g., the mean and covariance) of the system states, the 
computational cost can be significantly reduced. Moreover, in the EnKF there 
is no need to linearize the system model as in the EKF. Instead, the nonlinear 
problem becomes the one of how to approximate the statistics of the pdf of 
a Gaussian distribution which is transformed by a nonlinear function. This 
point of view leads to the idea of the unscented transform (UT) [HUT], as 
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will be introduced later. 



Depending on whether to perturb the observations or not, the EnKF can be 
classified into two types: stochastic and deterministic [T9|28] . The stochas- 
tic EnKF uses the observations and the corresponding covariance matrix to 
produce an ensemble of the disturbed observations, which is then used to up- 
date the background ensemble. For examples, see 0BEH1] • In contrast, the 
deterministic EnKF, often known as the ensemble square root filter (EnSRF 
hereafter), does not perturb the observations. Given a background ensemble, 
the EnSRF uses the observations to update the sample mean of the back- 
ground, while the analysis ensemble is produced based on the sample mean 
plus the perturbations. For examples, see [TfHISl] ; also see the review of [25] . 
Apart from the aforementioned EnKFs, there are also some other variants, 



In this paper we will introduce a framework of the EnKF incorporating the 
concept of the unscented transform [r5]TTT] , which will therefore be called the 
ensemble unscented KF (EnUKF for short). The basic idea of the EnUKF is 
that, at the end of the filtering step of an EnKF, one adopts the unscented 
transform to generate a set of carefully chosen system states, called the sigma 
points. The set of the sigma points is then treated as the analysis ensemble 
and propagated forward through the system model. At the next assimilation 
cycle, the mean and covariance of the background are estimated based on 
the propagations. In the Appendix, we show that, under the assumption of 
Gaussian errors, the accuracies of thus estimated mean and covariance are up 
to third order term in the Taylor series expansion (if available). Moreover, 
there are additional adjustable parameters in the framework which may be 
tuned to reduce the approximation error further. In contrast, for the ordinary 
EnKFs, we will show in the Appendix that the accuracies are normally up to 
second order. Therefore, by incorporating the unscented transform, one may 
improve the performance of the EnKF. 



This paper is organized as follows. In section [2] we will review the general 
framework of the EnKF. Moreover we will also introduce the idea of the un- 
scented transform. The accuracy analysis of the unscented transform is given 
in the Appendix. For comparison, we will also provide the accuracy analysis 
of the ordinary EnSRF. In section [31 based on the concept of the unscented 
transform, we will propose a modification scheme to the EnKF. In section [4j 
we will adopt the 40-dimensional model in [22] as the testbed to examine the 
performance of the EnUKF and compare it to the ordinary EnSRF. Finally 
we will conclude the work in section [SJ 
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2 Background 



2. 1 The framework of the EnKF 



The EnKF is a Monte Carlo implementation of the Kalman filter, which in- 
herits the framework of the KF by nature. For illustration, we consider the 
following scenario. Suppose that we have an m-dimensional discrete dynamical 
system 

x fe+1 = M k ,k+i (x fe ) + u fe , (1) 

where denotes the m-dimensional system state at the instant k, Ai k ,k+i is 
the transition operator mapping to Xfc +1 , and is the dynamical noise, 
which is assumed to be independent of the system state and observation noise 
(see below), with zero mean and covariance 

The above dynamical system is then measured by an observer 7i k such that 

y k = U k (x*.) + v fc , (2) 

where is the p-dimensional observation at instant k, and v& is the obser- 
vation noise, which is independent of the system state and dynamical noise, 
with zero mean and covariance R^. For convenience of discussion, we also as- 
sume that there is an n-member ensemble of the analysis X^_ x = {x^_ x i : 
i = 1, 2 • ■ ■ , n} available at the end of the (k — l)-th assimilation cycle. For 
description, we split the framework of the EnKF into the propagation (or 
prediction) and filtering steps. 



2.1.1 The propagation step 

We define the following set 

X£ = {A, ■■ = M k>k+1 (x£_ M ) , i = 1, 2, • ■ ■ , n) 

as the predicted background ensemble at instant k, which is the set of propa- 
gations of the analysis ensemble X^_ : at the previous assimilation cycle. The 
sample mean x| and covariance P b k can be evaluated according to the following 
unbiased estimators 0- 

b 1 71 
X fc = ~ E X fc,i ' 
i=l 

1 71 rp 

Pfc = — — r E Ki - *£) «, - *&) + Qk ■ (3b) 
1 In some works, e.g., [33], the authors may choose other estimators. 
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In practice, the approximation covariance P b k need not be calculated. Instead, 
it is customary to compute 

pJ* E - x 9 («* «*) - n k (x^)) T , 

P " h =—r E («* «*) " W* (*£)) («* (x^) - 7Y fe (**)) , 
n 1 i=i 

which avoids the problem of linearizing 7i k when it is nonlinear, although 
under such circumstance the ordinary KF algorithm becomes sub-optimal. 
The Kalman gain is obtained through 

K fc = P^(P^ + R fc )"\ (5) 
where is the covariance of the observational error at step k. 

2.1.2 The filtering step 

With additional information from the observations, one can update the back- 
ground ensemble according to a certain analysis scheme. For example, for a 
stochastic EnKF, the analysis ensemble X£ = jx^ : % = 1,2, • • • , n\ can be 
generated according to 

x^ = y* ti + K k (y fe , - H k (x^)) , for i = 1, 2, • • • , n, (6) 

where y k ,i is an observation sample generated by the normal distribution with 
mean y^ and covariance R^. Correspondingly, the sample mean and covariance 
of the analysis ensemble can be calculated according to 

1 n 

*2 = -E*m. ( 7a ) 

^ = rrTEfe-^)( x l-^) • ( 7b ) 
n 1 i=i 

For an EnSRF, for instance, the ensemble transform Kalman filter (ETKF) 
[1122], the analysis ensemble is generated by adding perturbations to the en- 
semble mean. On one hand, the ensemble mean x^ is obtained as follows 

n = ** + K k (y k - H (4)) . (8) 

On the other hand, let St be an m x n square root matrix of the background 
error covariance P| such Pj = S£ ( S >)B then an ,n x „ square root ma- 

2 In general can be obtained through some numerical algorithm, e.g., Cholesky 
decomposition. 
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trix of the analysis error covariance P% can be updated from S b k according 
to 

S a k = S b k T k , (9) 
where T k is an n x n transformation matrix derived in [I]. Thus the analysis 
error covariance P k is computed by 

P a k = s a k (s a k ) T . (10) 



Moreover, given x£ and S£, the analysis ensemble X£ = : i — 1, 2, ■ • ■ , nj 
is generated according to 

x^ = x£ + V^l (SI), , i = l,2, ■■■ ,n, (11) 

where (S^ denotes the i-th column of the square root matrix S£. After the 
analysis ensemble is generated, one propagates it forward to obtain the back- 
ground ensemble at the next instant and starts a new assimilation cycle. 



2.2 The unscented transform 



Given an m-dimensional Gaussian random variable x with mean x and covari- 
ance P xx , one problem of interest is how to estimate the mean and covariance 
of the transformed variable y = f (x), where f is a nonlinear vector function. 

In the ordinary EnKF, given a set of samples {xj,z = 1,2, ■ ■ ■ ,n} of the 
random variable x , the mean and covariance of the transformed variable y 
are estimated according to Eq. (jHJ), i.e., 

1 n 

y = -£ f (*0. ( 12a ) 

i=l 

1 n 

Pyy = -—7 E ( f (*) - y) ( f (^) - y) T • ( 12b ) 

n 1 i=X 



The unscented transform [16117] is a different method for the estimation 
problem. Suppose that x and P xx of the random vector x are unknown, 
but the sample mean x and covariance P xx can be estimated from the set 
{xj, i = 1,2, ■ ■ ■ , n}. Then a set of 2L + 1 state vectors i = 0, 1, • • • , 2L}, 
called the sigma points [HUT], can be generated according to 

X = x, 

Xi = x+(y/(L + \)P xx ^j , i = l,2,- ■■ ,L, ^ 
Xi = -k- (y/(L + \)P X3 ?) ^ , i = L + 1, L + 2, • • • , 2L, 
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where (^J (L + \)P XX ^J denotes the i-th column of the square root matrix 
{L + X)P XX , and A is an adjustable scaling parameter [TBUTj. 

Moreover, in the unscented transform, it also allocates a set of weights {Wi,i = 
0,1,..- ,2L} 

W 



L+ l (14) 

to the sigma points. In this way, the weighted mean and covariance of the 
discrete distribution {Xi, i = 0, 1, • • • , 2L}, 

2L 

x = ^w l x i , 

( 15 ) 

P xx = £ Wi (Xi - x) (xt - x) T , 

are exactly the same as x and P xx . If x follows a Gaussian distribution, then 
A can be chosen as A = 3 — L in order to let the first four weighted moments of 
the discrete distribution {Xi, % = 0, 1, • • • , 2L} match those of the distribution 
of x H5H7I. 



Because of the symmetry in the sigma points, the rank of the matrix Pxx is 
L. To avoid rank deficiency in the sample covariance matrix, it is suggested 
that the number of the sigma points be larger than twice the dimension of the 
vector x [MITT?] , or equivalently, L > m. 

After the transformation, the weighted sample mean and covariance of the 
transformed sigma points are calculated according to 

2L 

y = y £W i f{X i ), (16a) 

2L 

Pyy = £w« (f {Xt) - y) (f W - f) T + P (f (*o) - y) (f TO - y) T , 

j=0 

(16b) 



where the second term on the rhs of Eq. (116bj) is introduced to reduce the 



approximation error. In the case that x follows a Gaussian distribution, the 
choice of j3 = 2 is shown to be optimal [T7] . 

To take into account the effect of the model error in a more general situation, 
wherein the system model is described by y = f (x, u) (such that the model 
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error u is possibly not additive, but is assumed to follow a Gaussian process 
with mean zero and covariance Q), it is customary to adopt the joint state 
z = [x T , u T ] T and the system model is changed to y = f (z). In this case, the 
sigma points can be generated according to 









T 


Zi 


= 3> + (\ 




- x)P zz 


Zi 




/(L- 





1,2,--- ,L, 



(17) 



L + LL + 2,--- ,2L, 



where means the zero vector, and P zz is the sample covariance matrix of z 
such that 
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Q 



We will leave the awkward analysis of the estimation accuracies of both the 
ordinary EnKF and the unscented transform to the Appendix, where the gen- 
eral form y = f (z) described in terms of joint state is considered. We will 
show that, in the estimation scheme of the ordinary EnKF, the random sam- 
ples will generally introduce spurious modes in the transformed distribution 
even if the set of sample points has the correct mean and covariance [T6][TT] . 
while for the unscented transform, by carefully choosing the samples (i.e. the 
sigma points), the effect of sample error can be reduced. This fact leads to 
our argument that incorporating the unscented transform may benefit the 
performance of the EnKF. In the next two sections we will firstly introduce 
the modification scheme which incorporates the unscented transform for large- 
scale problems, and then proceed to conduct some numerical experiments to 
examine its performance. 



3 The ensemble unscented Kalman filter 



For large-scale problems, it is not practical to fulfil the requirement that 
the number of the sigma points should be larger than twice the degrees-of- 
freedom of the system model [12] . Therefore we will introduce some mod- 
ifications to make the unscented transform applicable for those problems. 
We will call the EnKF equipped with the modifications the ensemble un- 
scented Kalman filter (EnUKF). For consistency, we again take Eqs. (pP) and 
([2]) as the m-dimensional system model and the p-dimensional observer re- 
spectively. Moreover, we also assume that there is a set of sigma points 
{^k-n,i = 0, 1, • • • ,24-i} available at the (k — l)-th step, which is asso- 
ciated with the weights W^s given by Eq. (THj) . If there is only an "ordinary" 
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ensemble of the analysis available, one may first compute the sample covari- 
ance, and then conduct the truncated singular value decomposition (TSVD) 
[T5] to construct a set of sigma points, as to be explained later. 

With the zero mean of the dynamical noise, at the fc-th cycle, there ex- 
ists a set of the predictions of the propagated sigma points {X ki '■ X ki = 
M.k-i,k (^k-i,i) , * = 0, 1, • • • , 2£fc_ 1 } ) which is associated with a set of weights 

{w k 

-1,0, •"■ > Wjfc-i.aij.-i} obtained from the previous cycle (to be discussed 
later). The weighted sample mean and covariance are given by 

2//b-i 

*fc= £ Wik-i,iA* • (19a) 

i=0 

Pi = £ w^y (aj, - 4) {Ki " **f + (19b) 

i=0 

/3(4o-4) (x£ fi -*if + Qk. 



Moreover, the above evaluation scheme is also applied to the projection of the 
background ensemble such that 



2h 



y b k = £ w k ^n k (xl h 



i=0 



p xh = £ Wk-i,* - x») (H k (*£,,) - y b k ] 

,t (20) 



i=0 

+ (3 (x b k>0 - **) (H k (4 ) - y£ 

PL = £ w k ^ (n k (x b J - y b ) (n k (a* ) - f k 

i=0 

+ P (n k (4 ) - y b ) (n k (4 ) - y y T . 



For numerical reason, it is customary to re-write the above error covariances 
in terms of some square root matrices. To this end, we introduce two square 
roots, S£ and S£, which are defined by 



Sx 
k 



Sit 
k 



k-l,2h 



ir, 



fc-1,0 



yfc 



(21a) 



(21b) 



W k _ lt 2i k _ 1 [H k {X k!2lk l 



Yk 
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where W k _ 1Q = Wk-1,0 + P- Then the covariances can be re-written as 

P> = S£(S£) r + Q fc , (22a) 
P k xh = St(S k f, (22b) 
PL = S^(S^ T , (22c) 

Again, the Kalman gain follows Eq. (131), i.e., 

K fc = P^(p^ + R fc ) _1 . © 

With the above information, the mean and covariance of the analysis can be 
computed according to 

x£ = x 6 fc + K fc (y fc -?4(x£)), (23a) 
Pfc = P£-K fc (py T - (23b) 

Apart from obtaining the updated sample mean and covariance, we also aim 
to generate a set of sigma points as the analysis ensemble, which will then be 
propagated to the next assimilation cycle. For this purpose, one may consider 
using an existing EnKF scheme, for example, the ETKF. However, in order 
to avoid doubling the ensemble size at each assimilation cycle, some sigma 
points have to be discarded. To do this, the sample mean can be preserved 
by maintaining the symmetry about x£ among the remaining sigma points, 
while the corresponding sample covariance, denoted by P£, can only be an 
approximation to Pj£. This may appear to be a complicated problem for the 
existing EnKFs to design a selection criterion, because the perturbations pro- 
duced by them have no indications of the relative importance for covariance 
approximation. 

To tackle the above problem, the truncated singular value decomposition 
(TSVD) [13] is adopted in this work, which is similar to the idea of using 
singular vectors to produce initial perturbations for ensemble forecasting ([6], 
[T8| ch 6], [30], [31], and the references therein). Specifically, to produce the 
sigma points, a singular value decomposition (SVD) is firstly conducted on 
P k . Suppose that P% can be expressed as 

P a k = E k B K (E k ) T , (24) 

where Dk = diag(<j| u ■ ■ • , a\ m ) is a diagonal matrix consisting of the eigen- 
values cr^j's of P k , which are sorted in descending order, i.e., a ki > a k j > 
for i > j , and E K = [e k tl , • • • , e fc m ] is the matrix consisting of the correspond- 
ing eigenvectors e^'s. Next, a set of perturbations, generated in terms of the 
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first Ik values of cr^e^i, are added to the sample mean to form lk sigma 
points. Another Ik symmetric sigma points can be produced by subtracting 
the perturbations from the sample mean. Overall, in analogy to Eq. (ITS]) , the 
above procedure can be summarized as follows 

*k,i = + (h + A) 1 / 2 £7 fc)i e fc)i , z = 1, ■ • ■ , l k , (25) 
^fe.i = x-k ~ (h + A) 1 ^ 2 crfe i j_; fc efc i i_i fe , z = Zfe + 1, • • • , 2lk, 

where A is the adjustable scaling parameter. Using the generated sigma points 
as the analysis ensemble, the EnUKF is clearly an unbiased ensemble filter 

EDI. 



It is worthy to note that, Eq. (1251) does not require the full spectra of the 
eigenvalues and eigenvectors. Therefore, to reduce the computational cost in 
large scale problems, some fast SVD algorithms, e.g., the Lanczos or block 
Lanczos algorithm [TUJ ch 9], can be adopted to compute the first l k pairs of 
eigenvalues and eigenvectors (for example, see [29J). 

For convenience, we will hereafter call l k the truncation number. The choice of 
If. is important to the performance of the EnUKF, because it not only controls 
the number of sigma points to be produced, but also determines the quality 
of matrix approximation. Indeed, via SVD the matrices ~P% and can be 
expressed as 

m 

p* = E<< e *.<( e *.<) T . 

i=l 

respectively. It is clear that, if l k is too small, some important structures of P£, 
in terms of a\ fik^i ( e fc,«) T f° r i > hi w ^ be lost. However, as the computational 
cost is also a concern, it is not desirable for to get too large. Moreover, in 
many situations, if is large enough, a\ x may be already very small compared 
to the leading eigenvalues. Thus the improvement obtained by increasing If. 
becomes negligible. In this sense, one may choose a modest value of If. to 
achieve a tradeoff between accuracy and efficiency. 

In our implementation, we let lk be an integer such that 

k ■ 



a k i > trace (P£) /h k , i = 1, • • • , Z„ , 

\ (27) 
a\ i < trace (P a k ) fh k ,i>l k + l, 



where h k is the threshold at the k-th cycle (we will discuss how to choose h k 
in section H?Ti) . This is equivalent to saying that we construct the sigma points 
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based on the eigenvectors such that their corresponding eigenvalues are larger 
than a specified tolerance. Moreover, to prevent l k getting too large or too 
small, we also specify a lower bound li and an upper bound l u and adjust the 
threshold h k so that li < l k < l u . 

Under the assumption of Gaussian error distribution, it can be verified that the 
perturbations of the sigma points, in terms of (ifc+A) 1 ' 2 ^^^ for i — 1, ■ • • , l k , 
are equally likely in the sense that their values of the probability density 
function (PDF) 

p(*x) = (27r) m/2 (detP»)-^exp |-i (5x) T (p")" 1 (5x)} (28) 

are the same (also see the discussions in [33]), where det • means the deter- 
minant of a matrix. Therefore it is natural to assign an identical weight to all 
the perturbations. Consequently, in the spirit of Eq. (JHJ), a set of weights can 
be constructed as follows 



h + y ( 29 ) 

W k,i = 777] , \T> * = 1> ' ■ • , 2l k- 
I [L k + A) 

Finally, all the sigma points in Eq. (1251) . which are associated with a set of 
weights given in Eq. (1291) . are propagated forward to the next assimilation 
cycle. 

We adopt the time averaged relative rms error (relative rmse for short) to 
measure the performance of the EnUKF, which is defined as 



A- 



max 



J2\\±t- x th/\Kh, (so) 



where k max is the maximum assimilation cycle, x^T denotes the truth (the state 
of a control run) at the k-th cycle, and ||«||2 means the L2 norm. 

Moreover, we also use the time averaged rms ratio to examine the similarity of 
the truth to the sigma points, which also qualitatively reflects the performance 
in estimating the error covariance, e.g., overestimation or underestimation (cf. 
[Tf3~4] and the references therein). As an estimation, the time averaged rms 
ratio, denoted by R, is computed by [Tf3l] 



R=T^-J2 Wk + i)\\n-^h/EK,i-**h, (3i) 

Kmax k= i i=0 

while the expectation of the rms ratio is [Tf3l] 

Re = yfiUff + l)/(2Ze// + 1) , 
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where l e ff is the "effective" truncation number over the whole assimilation 
window. Hence, if the truth is statistically indistinguishable from the sigma 
points, the values of R and R e shall be very close. Note that R e ~ 0.71 for 
any large l e ff, so for simplicity we let l e ff equal the average of the truncation 
number I, i.e., i e // = I — J2i=± x lk/k max . R > R e means that the covariance 
of the sigma points underestimates the error of the state estimation, while 
R < R e implies the opposite, i.e., overestimation of the error of of the state 
estimation [24T34] . 



4 Numerical experiments with a 40-dimensional system 

This section is dedicated to examining the performance of the EnUKF through 
the numerical simulations, and studying the effects of the parameters on the 
performance of the EnUKF. To this end, we choose the m-dimensional sys- 
tem model due to Lorenz and Emanuel [2~T|f2"2"] (LE98 model hereafter) as the 
testbed. The LE98 model is a simplified system proposed to model atmo- 
spheric dynamics, which "shares certain properties with many atmospheric 
models" [22]. We consider the perfect model scenario, wherein the governing 
equations are described as follows 

^ = (x i+1 - Xi_ 2 ) Xi-i - %i + F, i = 1, • • • , m. (32) 
at 

The quadratic terms simulate the advection, the linear term represents the 
internal dissipation, while the constant F acts as the external forcing (|21j). 
The variables a^'s are defined cyclically such that X-i = x m _i, xo = x m , and 

Xm+l X\. 

We choose the observer Tik to be a time-invariant identity operator. Specifi- 
cally, given a system state x^ = [xk,i, • • • , Xk, m ] T at the fc-th assimilation cycle, 
the observations are obtained according to 

y fc = H k (x k ) + v fc = x fc + v fc , (33) 

where follows the m-dimensional Gaussian distribution N(0, R&) with the 
covariance matrix being the m x m identity matrix I m . 

Note that in Eq. fl30l the measure e r can also be interpreted as the time- 
averaged noise level of the trajectory {x£}^™^ (with respect to the truth). 
From this point of view, we can define the divergence of a filter in a restric- 
tive sense: Suppose that the relative rmse of the observations is e° bv , which, 
with the identity observation operator TCk, is defined as 1 1 y ^ — x|T 1 1 2 / 1 1 x | r 1 1 2 = 
||vfc|| 2 /||x| r || 2 . If e r > e° bv , then we say the filter is divergent because in such 
circumstances, the trajectory {x^}^™^ obtained by the filter, on average, is 
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more noisy than the observations, which implies that it might not make any 
sense to use the filter for assimilation. 

In our experiments, we set m = 40 and F = 8 and integrate the system 
through the fourth-order Runge-Kutta method. We choose the length of the 
integration window to be 100 dimensionless units, and the integration time 
step to be 0.05 units (corresponding to about a 6-h interval in reality, see 
[22]). thus there are 2000 assimilation cycles overall. 

For the LE98 model, in the numerical experiments (not reported there) we 
found that, the EnSRF, implemented following either work of [H4l34] . out- 
performs the stochastic EnKF, which is consistent with the result reported 
in [33], while the performances of the EnSRFs are very close. Therefore, here 
we choose to compare the EnUKF with the EnSRF only. Moreover, for an 
ordinary EnSRF, we will introduce the spherical simplex centering scheme 
to its analysis ensembles. The reason to do this is two-fold. Firstly, with the 
spherical simplex centering scheme, the produced analysis perturbations are 
equally likely in probability under the assumption of Gaussian error distri- 
bution [33]. Secondly, with the centering scheme, an EnSRF can be shown 
to be an unbiased ensemble filter, which avoids an error that systematically 
underestimates the analysis covariance (see [20] for the details). In this work, 
we single out the ETKF [1] for our experiments, and equip it with the spher- 
ical simplex centering scheme following the work [33], which will be further 
discussed below. 



4-1 Issues in implementing the EnUKF and the ETKF 

When implementing the EnUKF, an issue worth of special attention is the 
positive semi-definiteness of the covariance matrices. Normally, we require 
Ik + A > so that in Eq. fT29l . the weights Wk/s are positive for i > 0. 
Nevertheless, the weight Wk,o can be negative if A < 0. If it is so, when 
computing the background covariances according to Eqs. (I19bl) and (1201) . the 
positive semi-definiteness may not be guaranteed. However, one may note 
that, in Eqs. (I19bl) and (I2U1) . the effective weight of X^ +1 is actually Wk,o + ft 
{ft > 0). So in order to guarantee the positive semi-definiteness, we shall 
choose the parameters A and ft properly to satisfy that Wk,o + ft > and 
Ik + A > 0. Given Wk,o = A/ (Ik + A), it means that A > —fth/ (1 + ft)- Since Ik 
is bounded within [//, l u ], by letting A > —fth/ (1 + ft), one can guarantee the 
positive semi-definiteness. 

The threshold hk in Eq. (127|) is chosen in the following way. At the beginning we 
specify a threshold h\. If hi is a proper value such that l\ satisfies \\ <l\ < l u , 
then we keep hi and at the next cycle we start with h 2 = h±. If hi is too 
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small such that l\ < li, then we replace hi by l.lhi + 200. We continue 
the replacement until l\ falls into the specified range, or the number of the 
replacement operations is up to 30 (in this case we simply put l± = //, regardless 
of what hi is). Similarly, if hi is too large such that li > l u , then we replace 
hi by h\/l.l — 200, we continue the replacement until li falls in the specified 
range, or the number of the operations is up to 30 (in this case we simply put 
li = l u ). After the adjustment, at the next cycle we start with h 2 = hi and 
adjust it (if necessary) to let l 2 fall into the specified range, and so on. In this 
way, one can obtain the threshold at each cycle. 

To apply the spherical simplex centering scheme to the ETKF, we follow the 
algorithm proposed in [17] to construct a centering matrix U, where U follows 
Eq. (C15) in [33]. Compared to the alternative centering matrix provided in 
[33] , we favor the one described in Eq. (C15) because it is time- invariant. More- 
over, its construction is independent of the concrete system in assimilation, 
and does not involve the operation of matrix inversion and the observation 
operator. 

In our experiments, we process the observations simultaneously in both the 
EnUKF and the ETKF. In order to improve the performances of the filters, 
we also consider two additional techniques. One is the method of covariance 
inflation, which is based on the observation that the covariance of the analysis 
error will be systematically underestimated in the EnKF [M]. Therefore, it 
can be beneficial to increase either the background error covariance, or the 
analysis error covariance [2"f 251134] . In this work, we follow the method used 
in [2f34"] and choose to multiply the perturbations to the sample mean x£ of 
the analysis by a constant 1 + 5, which is equivalent to increasing the analysis 
error covariance by a factor (1 + 5) 2 . Note that, if one chooses to process the 
observations in a serial way (e.g. [31]) with a covariance factor 5 S , then after 
processing an m-dimensional observation, the corresponding covariance will 
be increased by a factor of (1 + 5 s ) 2m . The relationship between the inflation 
factors 5 and 5 S is thus given by 5 = (1 + 5 s ) m — 1, therefore one will find that 
the inflation factor adopted in this work is much larger than those used in, for 
example, [3"4"j . 

The other technique is covariance filter [TTlfTS] . which introduces the Schur- 
product to a covariance matrix in order to reduce the effect of sample errors. 
We say a length scale of covariance filter is optimal within a certain range if 
it minimizes the relative rmse among the values in test. As an example, in 
Fig. [T] we plot the relative rms errors of the EnUKF (upper panel) and the 
ETKF (lower panel) vs the length scale. Both filters start with the same initial 
condition, and takes no covariance inflation (i.e., 5 = 0). The length scale is 
varied from 40 to 400, with an even increment of 40 at each step. The other 
settings of the filters are as follows: For the EnUKF (upper panel), the initial 
ensemble size is 4. (3 = 2, A = —2, the lower bound k = 3, the upper bound 
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l u = 6, and the threshold h\ = 1000; For the ETKF, the ensemble size is 13. 
From Fig. [IJ it is clear that the optimal length scale of the EnUKF is 200, 
while the optimal length scale of the ETKF is 240. For simplicity, we choose 
l c = 240 for both filters in the subsequent simulations. 



4-2 Comparison between the EnUKF and the ETKF 



For comparison, we randomly select an initial condition xi, and use it to start 
a control run. The observations Y = {yk}t=i are obtained by adding Gaus- 
sian white noise to the states of the control run at each cycle, in accordance 
with Eq. (I33l . In subsequent simulations, both the EnUKF and the ETKF will 
start with the same initial condition xi, and use the same observations Y for 
assimilation. To initialize the filters, at the first assimilation cycle we randomly 
generate a background ensemble = {x^j : i — 1, • ■ ■ , n}. Given and x 1; 
the ETKF is already able to start running recursively. For the EnUKF, how- 
ever, at the first cycle there is no propagated sigma points from the previous 
cycle. So, similarly to the ETKF, one may use the background ensemble Xj to 
compute the sample mean and covariance of the analysis, and then generate 
the sigma points accordingly. After propagating the sigma points forward, the 
EnUKF can start running recursively from the second cycle. 



For the EnUKF, we let the parameters (3 = 2, A = —2, the threshold hi = 
1000, the lower bound l t = 3, the upper bound l u = 6, the length scale of 
covariance filter l c = 240, and the covariance inflation factor 8 vary from 
to 10, with an even increment of 0.5. We consider the scenarios with different 
ensemble sizes n = 3,4, 5,6 at the first assimilation cycle in order to explore 
the effect of initial ensemble size on the performance of the EnUKF. The 
corresponding relative rms errors and ratios, as functions of the covariance 



inflation factor, are plotted in Figs. 2(a) and 2(b) respectively 



From the above two figures, it can be seen that different initial ensemble sizes 
n = 3, 4, 5, 6 leads to similar behaviors of both the relative rmse and rms ratio. 
Interestingly, a larger initial ensemble size does not necessarily guarantee a 



smaller rmse error. This can be observed either from Fig. 2(a) by fixing the 
covariance inflation factor 8 at some point, say 8 = 1.5, or from Table [TJ by 
comparing the minimum relative rms errors. 



Fig. 2(a) shows that, as the covariance inflation factor 8 increases from 0, the 
relative rmse of the EnUKF tends to decline. However, if 8 gets too large, say 
8 > 6, then further increments in 8 will instead boost the relative rms errors. 
An examination on the rms ratio also reveals the same trend, although, as 
indicated in Fig. 2(b), the turning points, now at 8 = 7.5 for n = 3,4, 6 and 



8 = 7 for n = 5, are larger than those of the relative rms errors. To make 
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the sigma point indistinguishable from the truth (i.e., rms ratio ~ 0.71 ), one 
needs the inflation factor 5 ~ 2.5. However, modestly larger inflation factors, 
say, 2.5 < 5 < 6, can benefit the performance of the EnUKF in terms of 
the relative rmse, although they also cause the over-estimations of the error 
covariances. 



For the ETKF, we let the length scale l c of covariance filter and the covariance 
inflation factor 5 be the same as those in the EnUKF. Suppose that in a run of 
the EnUKF we have the average truncation number /. Then for comparison, 
we consider the ETKF with an ensemble size n = ceil (21 + 1), where ceil(s) 
means the nearest integer that is larger than, or equal to, the real number s. In 
our experiments, the EnUKF with different initial ensemble sizes n = 3, 4, 5, 6 
leads to the same value ceil(2Z + 1) = 13, so it is not surprising to find in 



Figs. 3(a) and 3(b) that the relative rms errors and ratios of the ETKF, 
which correspond to the EnUKF starting with different initial ensemble sizes, 
actually coincide. 



From Fig. 3(a), one can see that, starting from § = 0, as the covariance 
inflation factor increases, the relative rmse of the ETKF tends to decrease. 
However, unlike the situation in the EnUKF, in the test range, as 5 gets larger, 
say 5 > 7, the corresponding relative rmse enters a plateau region. The rms 
ratio indicates a similar behaviour. As 5 increases, the change of the rms ratio 
becomes smaller, or in other words, the curve appears more and more flat. In 
order to make the ensemble in the ETKF indistinguishable from the truth, one 
needs the inflation factor 5 ~ 1.5. Like the EnUKF, modest overestimation 
of the error covariance (i.e., 5 > 1.5) can also benefit the performance of the 
ETKF in terms of the relative error. 



We use the minimum relative rms errors of the EnUKF and the ETKF to 
compare their performances. To this end, in Table [1] we list the minimum 
relative rms errors of the EnUKF with different initial ensemble sizes, and the 
minimum relative rms errors of the ETKF with the ensemble size about twice 
the average truncation number plus Moreover, reference, we also show 
the average noise level (relative rmse) in the observations. 

From Table [H one can see that the average noise level in the trajectory consist- 
ing of the observations is roughly 0.226 (22.6%), while the minimum relative 
rmse (i.e. noise level) of the assimilated trajectory through the ETKF scheme 
is about 0.207 (20.7%), a reduction of 1.9% noise level compared to the ob- 
servations. Similarly, the minimum relative rmse of the assimilated trajectory 
through the EnUKF scheme is roughly less than 0.175 (17.5%), a reduction 
of more than 5% noise level compared to the observations, and more than 3% 

3 Note that different initial ensemble sizes n = 3, 4, 5, 6 in the EnUKF lead to the 
same ensemble size in the ETKF. Therefore, in Table [TJ the ETKF has the same 
minimum relative rmse in different rows. 
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compared to the ETKF scheme. In this sense, both the ETKF and the EnUKF 
do not diverge as their minimum relative rms errors are less than the average 
noise level in the observations, yet the EnUKF exhibits a better performance 
in state estimation than the ETKF. 



4-3 Effects of parameters on the performance of the EnUKF 

There is a set of adjustable parameters in the EnUKF, e.g., threshold hi, 
lower bound li and upper bound l u , A in Eqs. f l25l) and fl29l) . and j3 in Eqs. 
( 119bl) and (|20j) . In this section we study the effects of these parameters on 
the performance of the EnUKF. Note that in the previous section, we have 
already examined the effect of the inflation factor on the performance of the 
EnUKF. Thus in the subsequent experiments, we will just fix the inflation 
factor at a particular value (say, zero, but other choice will also do) for the 
sake of simplicity. 

4-3.1 Effects of the threshold hi and the bounds l\, l u 

The parameters hi, li and l u determine the truncation numbers l^s. To ex- 
amine their effects on the performance of the EnUKF, we fix the covariance 
inflation factor 5 = 0, the length scale l c = 240, parameters A = —2 and (3 = 2. 
We specify the upper bound l u = 6 in the experiments, but vary the lower 
bound li from 3 to 6. This choice is used to represent the typical situation in 
data assimilation, wherein the ensemble member is often much lower than the 
dimension of the dynamical system. We also vary the threshold hi, in the log- 
arithmic scale log 10 hi, from 2 to 5.5, with an even increment of 0.5 each time. 
This range represents the moderate values of hi so as to make the truncation 
numbers neither too large nor too small. In all the experiments, we initialize 
the EnUKF with the same initial conditions and background ensemble, with 
the ensemble size n = 4 at the first assimilation cycle. 

In Fig. H] we show the simulation results. Clearly, the larger the threshold hi 
and the bound l\ are, the larger the truncation number l\~ tends to be, which, 
however, does not mean the better performance in terms of the relative rmse. 
Indeed, in Fig. HI there exists the same optimal threshold log 10 hi = 3 for 
lower bounds l\ = 3,4,5 0, while the thresholds larger than this value result 
in larger relative rms errors. For the lower bound li = 6, its relative rms errors 
are smaller than, or at least equal to those of the bounds k = 3,4,5 in most 
cases. However, at log 10 hi = 3, the relative rmse given li = 6 is worse than the 



4 h = lu = 6 means I}. = 6 at every cycle, while the threshold hi does not affect 
the choice of I},. 
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other cases. To explain these phenomena, we conjecture that, a too small trun- 
cation number Ij, is not likely to achieve a performance as good as a modest 
value because it means poor quality of covariance approximation. In contrast, 
a too large truncation number also does not necessarily achieve a better per- 
formance than a modest value because, at some local points of the attractor, 
a too large truncation number may introduce some spurious structures-from 
the null space of the SVD- into the sigma points, which are then treated as 
equally likely as the other sigma points, and propagated forward to the next 
cycle. The effect of the spurious structures can be accumulated and eventually 
deteriorates the overall performance. 



4-3.2 Effects of the parameters A and (3 



We proceed to examine the effects of the parameters A and (3. In the experi- 
ments, we set the covariance inflation factor 5 = 0, the length scale l c = 240, 
lower bound li = 3, upper bound l u = 6, threshold hi = 1000, and initial en- 
semble size n = 4. We consider four scenarios with (3 = 0,2,4,6 respectively. 
For each value of (3, we compute 20 values of A. To guarantee the positive 
semi-defmiteness of the sample covariances, we start A from —f3li/(l + j3), 
with an even increment AA = 1 each time. In particular, when j3 = and 
A = 0, the effective weight W^q + f3 of the ensemble mean x£ equals zero for 
any k. Therefore, in this case, the unscented transform is actually equivalent 
to the analysis scheme of positive-negative pairs (PNP) in the literature (cf 
[33] and the references therein). 

We plot the simulation results in Fig. As can be seen, when (3 increases from 
to 6, the minimum relative rmse for a given value of (3 intends to decrease. 
This may be interpreted as follows: In Eqs. f|19bj) and (T201 . the second terms 
on the rhs also act like a covariance inflation technique. Therefore, similar to 
the covariance inflation factor 5, a larger value of j3 tends to result in a smaller 
minimum relative rmse. 



However, for each fixed f3, there is no clear trend of the optimal value of 
the parameter A. A larger value of A does not imply a smaller relative rmse, 
or vice verse. Particularly, from Fig. 5(a) it can be seen that, by choosing 



suitable values for the parameters A and (3, the EnUKF can outperform the 
EnKF equipped with the analysis scheme of positive-negative pairs. 



As an explanation of the above phenomenon, one may note that, with the 
other parameters fixed, A is the parameter that determines the relative weights 
between the sample mean and the other sigma points (cf. Eqs. ff25l) and (129]) ). 
If the underlying system state is linear, then in principle one shall be able 
to compute the optimal relative weights between the sample mean and the 
other sigma points under the Gaussianity assumption, and thus determine 
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the optimal value of A. Nevertheless, the existence of nonlinearity may make 
the problem intractable. For nonlinear systems, the optimal relative weights 
(hence the parameter A) may vary from cycle to cycle. However, to search for 
the optimal parameter A at each assimilation cycle will be computationally 
expensive, thus in our ex per iments, we chose to fix the parameter A in the 
same assimilation window 5 1. This choice cannot reflect the variation of the 
optimal values of A at different assimilation cycles, therefore it would not be 
surprising to see that there is a lack of clear trend of the optimal value of A 
in Fig. (5(a) ). 



5 Conclusions 



A new ensemble Kalman filter scheme, called the EnUKF, was introduced 
in this work. The ensemble generation scheme of the EnUKF is similar to 
the idea of positive-negative pairs (PNP) ([33] and the references therein), 
but it differs from the PNP scheme in that, apart from generating symmet- 
ric positive-negative pairs, the EnUKF also propagates the ensemble mean 
forward. Like the EnKF equipped with the PNP scheme, in the EnUKF all 
symmetric perturbations are associated with an identical weight. Nevertheless, 
a different weight can be assigned to the ensemble mean. In this sense, the 
EnUKF can be deemed as a hybrid of central forecast (i.e., the forecast made 
by propagating the ensemble mean solely) and ensemble forecast (i.e., the 
forecast made by propagating the ensemble forward), while the PNP scheme 
is a special case of the unscented transform with the weight of the ensemble 
mean being zero. 

From the analytic results in the Appendix, it can be seen that, in estimation 
of the sample mean and the covariance of a transformed random variable, 
the accuracies of the EnUKF will be at least up to second order, no matter 
whether the original random variable follows a Gaussian distribution or not. 
If the original random variable follows a symmetric distribution (not necessar- 
ily Gaussian), then the accuracies of the estimations increase to third order. 
Moreover, additional parameters, such as A and f3, are available in the un- 
scented transform to improve the estimation accuracies by choosing proper 
values. 

By comparing the Taylor series of the transformation function term by term 
in the Appendix, we also show that the EnUKF has better accuracies than the 
ordinary EnKF. For numerical verification, using the LE98 model, we com- 
pared the performances of the EnUKF and the ETKF in terms of the relative 

5 Here by "assimilation window" we mean the time window from the first assimi- 
lation cycle to the maximum. 
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rms errors. Experiment results confirmed that incorporating the unscented 
transform into an EnKF can benefit its performance. 
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Appendix: Accuracies of the sample mean and covariance of the 
EnKF and the unscented transform 



In this Appendix we analyze the accuracies of the ordinary EnKF and the 
EnUKF in estimating the mean and covariance of a random variable trans- 
formed by a nonlinear function. In analysis, we assume that the nonlinear 
function can be expanded in a Taylor series, which converges to the true value 
of the transformation [TdITTT] . 



The actual mean and covariance of the transformed variable in terms of Taylor 
series 



Given a vector z and a Gaussian perturbation 5z with zero mean and covari- 
ance P 22 , let us first expand the transform y = f(z + Sz) in a Taylor series 
around the point z, then the mean y is given by 



y = E(f(z + <5z)) 




(34) 



where the operator 



D 5z = dz 1 V. 



(35) 



Similarly, the covariance matrix P yy is given by 



P 



.y.v 



E ((y - y) (y - y) T ) 



jp zz r + e 




3! 



DLf (DLf f , Df z f (D 5z f) 



(36) 



+ 



2! x 2! 3! 




where J = (Vf) T 



z 



is the Jacobian matrix of f at z. 
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Accuracies of the EnKF 



In the EnKF, given an ensemble {zj}™ =1 , the sample mean, in terms of Taylor 
series, is given by 



. J27=i d| z i 127=i f 

y =f (z) + + 



+ 



n x 2! 
ELiDLf 



nx3! 



(37) 



+ 



n x 4! 

where = z$ — z with z = X)^=i z i,M being the sample mean of z. 
In the rhs of Eq. ([3] 



sr=lI V = iv T f-E^zf) vf, 

nx2! 2! V n & / 

which is a biased estimation of the second term on the rhs of Eq. (134)) . More- 
over, spurious modes may rise in higher order terms of Eq. fl37j) . for example, 
the third order term 



1 EoL l f = v T (l-tszMvszj) Vf 



72 , 
«=1 



(3? 



in general will not vanish. 
Similarly, we have the sample covariance 

1 



P -TP TU 

yy — zz° 



+ 



[n - 1) x 2! ^ 



E D fal f(D? Bi f)- 



n — 1 



V 



3! 

r 



^i% f2 K, f ) , Er =1 DL,f(p^f) G 



(39) 



2! x 2! 



n — 1 



V r P, 2 V 
2! 



3! 



2! 



Note that here 



u 1 i=l 



Comparing Eq. (1391) and Eq. (1361) . we note that 
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Compared to Eq. (1361) . there are also spurious modes arising in higher order 
terms in Eq. (1391) ; 



There is a bias in the term 



V T P 22 V f 



V J P 22 V f 



Accuracies of the unscented transform 



For the unscented transform, given a set of sigma points {Zi}fto w ith mean 
z and covariance P zz , the sample mean, in terms of Taylor series, is given by 



2L 
i=0 



(40) 



Note that in Eq. (|40j) . the third order terms vanish because of the symmetry 
in the sigma points. 



Comparing Eq. (jjQl) with Eq. (1341) . we note that, second and third order terms 
in Eq. fj4*0l) match those in Eq. fl34l) exactly, while the difference starts from 
fourth order terms. 



Similarly, the sample covariance is given by 



2L 



mi 



EWiCf^-y) (f(^)-y) 



i=0 



z3PzzJT+ 2(L + X) 



E^D^f(DL,f) : 
3! 



Eiii D,5 Zl f 2 (p] z i) T t EL L iDLT(D 5zi f)^ 



2! x 2! 



'V T P ZZ V 



V T P„V 
2! 



3! 

l T 



(41) 



Clearly, unlike Eq. (1391) . there is no spurious mode rising in third order terms 
in the unscented transform. Moreover, there is also no bias in the term 



(V T P Z2 V) f 



V T P 22 V)f 



T 



To further reduce the approximation errors in fourth order terms, one may 
introduce an additional term, (3 (f (Z ) — y) (f (Z ) — y) T , to the sample co- 
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variance such that it can be re-written as 

2L 

Pyy = E Wi (f m - y) (f (Z,) - y) T + (3 (f (Z ) - y) (f (Z ) - yf . 

i=0 

dEbj) 

It can be shown that (cf. Eq. (21) of [H]), i n the case that z follows a Gaussian 
distribution, choosing (3 = 2 will minimize the approximation errors in fourth 
order terms. 
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Table 1 

Minima of the relative rms errors in Figs. 2(a) and 3(a) and the average relative 
rmse (noise level) e° bv of the observations. 



Initial ensemble size EnUKF ETKF Observations 



n 
n 
n 



3 
4 
5 



n = 6 



0.1719 0.2074 

0.1722 0.2074 

0.1730 0.2074 

0.1753 0.2074 



0.2256 
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Length scale of covariance filter 
(a) Relative rmse vs the length scale of covariance filter for the EnUKF 
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(b) Relative rmse vs the length scale of covariance filter for the ETKF 

Fig. 1. Effects of the length scale of covariance filter on the performances of EnUKF 
and the ETKF. From the figures, the optimal length scale of the EnUKF is 200, 
while the optimal length scale of the ETKF is 240. 
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(a) Relative raise of the EnUKF vs the covariance inflation factor 5 
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(b) RMS ratio of the EnUKF vs the covariance inflation factor 5 



2. Effects of the covariance inflation factor 5 on the performance of the EnUKF. 
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-o- Corresponding to the EnUKF with initial ensemble size = 3 
-> Corresponding to the EnUKF with initial ensemble size = 4 
-*- Corresponding to the EnUKF with initial ensemble size = 5 
-a- Corresponding to the EnUKF with initial ensemble size = 6 
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Covariance inflation factor 5 
(a) Relative rmse of the ETKF with the ensemble size equal to ceil (21 + 1) 
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-> Corresponding to the EnUKF with initial ensemble size = 4 
-*- Corresponding to the EnUKF with initial ensemble size = 5 
-a- Corresponding to the EnUKF with initial ensemble size = 6 
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(b) RMS ratio of the ETKF with the ensemble size equal to ceil(2/ + 1) = 13 



Fig. 3. Effects of the covariance inflation factor 5 on the performance of the ETKF. 
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Fig. 4. Relative raise of the EnUKF vs the threshold h\ (in the scale of log 10 ) with 
different lower bounds l\. 
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(d) Relative rmse vs A with (3 = 6 



Fig. 5. Effects of the parameters (3 and A on the performance of the EnlJKF. 
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